Nonextensivity of the cyclic Lattice Lotka Volterra model 
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We numerically show that the Lattice Lotka- Volterra 
model, when realized on a square lattice support, gives rise to 
a finite production, per unit time, of the nonextensive entropy 



! (Si = —Jjj Pi In This finiteness only oc- 
0.5 for the d = 2 growth mode (growing droplet), 



curs for q - 

and for q = for the d = 1 one (growing stripe). This strong 
evidence of nonextensivity is consistent with the spontaneous 
emergence of local domains of identical particles with fractal 
boundaries and competing interactions. Such direct evidence 
is for the first time exhibited for a many-body system which, 
at the mean field level, is conservative. 
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Many natural and artificial systems are known today to 
be hardly, or not at all, tractable within Boltzmann-Gibbs 
statistical mechanics, hence the usual thermodynamics. Such 
is the case of systems which include long-range interactions 
or long-range microscopic (or mesoscopic) memory, or other 
sources of (multi)fractality. Phenomena where many spa- 
tial and/or temporal scales are involved typically exhibit 
power laws. The celebrated Boltzmann-Gibbs (BG) entropy 
Sbg = — Pi In Pi appears to be inadequate for handling 
the thermostatistics associated with such situations. This is 
due to the fact that the corresponding stationary states do 
not emerge through ergodic dynamics. An ubiquitous class of 
the above anomalous systems have nonlinear dynamics which 
generate weak chaos, in the sense that the sensitivity to the 
initial conditions is less-than-exponential in time. Such sit- 
uations quite naturally accomodate with an entropy which 
generalizes the BG one, namely 



S q = 



EiP? 



(q £ ft; Si = Sbg)- 



G) 



For independent systems A and B (i.e., such that pf 3 +B = 
pfpf), this entropy satisfies S q (A + B) = S q (A) + S q {B) + 
(1 — q)Sq(A)S q (B). It is due to this property of nonexten- 
sivity that the thermostatistical formalism based on Eq. (1) 
is usually referred to as nonextensive statistical mechanics [1] 
(see [2] for recent reviews). This theory has received many 
applications in areas such as self-gravitating polytropes [3], 
electron-positron annihilation [4], turbulence [5], motion of 



Hydra viridissima [6], anomalous diffusion [7,8], classical [9] 
and quantum [10] chaos, long-range-interacting many-body 
Hamiltonians [11], option pricing [12], particular biological 
processes involving a large number of degrees of freedom 
[13,14], among others. 

The lattice Lotka- Volterra (LLV) model [15] is known to 
well mimic simple chemical reactions, predator-prey systems, 
and other biological and echological phenomena. It has re- 
cently been studied with cyclic interactions amongst three 
species, as a modification to the original Lotka- Volterra model 
[16]. The LLV has exhibited spatial clustering and fractal- 
ity, even when implemented on low dimensional lattices [17]. 
Since fractality is a distinct sign of nonextensivity, it is natural 
to pose the question whether the LLV model is indeed consis- 
tent with the nonextensive premises. The aim of the present 
paper is to verify that this model exhibits a strong and di- 
rect evidence of having g / 1. In particular we study the 
LLV model with various initial conditions, namely at a) the 
domain formation mode, b) the nucleus (or droplet) growth 
mode, c) the stripe growth mode and d) the roughening mode. 
Modes (b) and (c) enable, as we shall see, the direct calcu- 
lation of q. Several analytical or numerical calculations of q 
exist already in the literature, but this is the first time such 
evidence is directly provided, through the time evolution of 
S q itself, on a many-body system which is conservative at the 
Mean Field (MF) level. 

The LLV model is a minimal complexity model, with MF 
conservative dynamics which can be directly implemented on 
lattice and involves only two reactive species Xi and X2 (ad- 
sorbed on a lattice support) and the empty sites of the support 
S. All reactive steps are bimolecular and the reaction occurs 
via hard core interactions. Schematically, the LLV model has 
the following form [15]: 



Xi + x 2 
x 2 + s 

S + X x 



fcl 
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(2a) 
(2b) 
(2c) 
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In particular, a particle X\ adsorbed on a lattice site changes 
its state into X2 when it is found in the neighborhood of an- 
other X2 particle. This step (2a) is an autocatalytic reactive 
step. A particle X2 desorbs leaving an empty site S, if in the 
neighborhood another empty site 5* is found. This step (2b) 
is a cooperative desorption step. Finally, a particle X\ can 
be adsorbed on an empty lattice site S if in the neighborhood 
another X\ particle is found. This step (2c) is a cooperative 
adsorption step. 

We now recall briefly some of the mean field (MF) and lat- 
tice properties of the LLV, which have been studied in detail 
in previous works [15,17]. 
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In the MF approximation the LLV model, Eqs. (1), can be 
described via the kinetic, rate equations: 



= xi(-k s x 2 + k 2 s) 

at 

= x 2 (k s xi - kis) 

ds 

— = s(-k 2 xi + kix 2 ) 
at 



(3a) 
(3b) 
(3c) 



where Xi, x 2 and s correspond to the mean coverage of the lat- 
tice with particles X\, X 2 and empty sites S respectively. In 
Eqs. (3), the mean coverages satisfy identically the conserva- 
tion condition x\ + x 2 + s — C, where C is a constant which 
can be chosen equal to unity, corresponding to interpreting 
xi, x 2 and s as fractions of the overall lattice respectively oc- 
cupied by X\ particles, X 2 particles or being empty. Using 
C = 1 it is possible to eliminate one of the three variables, say 
s = 1 — xi — x 2 , and, to reduce system (3) to two equations. 
This reduced system admits four steady state solutions, three 
of which are trivial, and one non-trivial [15]: 



xis = 0, x 2s — (empty lattice) 

xis = 1, x 2s — (lattice poisoned by Xi) 

xis = 0, x 2s — 1 (lattice poisoned by X2) 
ki k 2 



Xls = 



ki + k 2 + k s 



x 2s 



ki + k 2 + k s 



(4a) 
(4b) 
(4c) 

(4d) 



A linear stability analysis shows that the trivial states are 
saddle points while the nontrivial one is a center compati- 
ble with an additional constant of motion C' = x^x^il — 
xi — x 2 ) ks [18] at the MF level. Fig. la depicts the temporal 
evolution of the system for typical values for (ki,k 2 ,k s ) and 
initial conditions. The black solid line represents the concen- 
tration of Xi and the dashed line the concentration of X 2 . 
The motion is periodic but non-harmonic. The amplitude 
of the periodic motion, for given parameter values, depends 
solely on the initial conditions [15,17]. At this level of descrip- 
tion the system size does not enter into the calculations since 
the MF approximation involves only average concentrations. 

To mesoscopically describe the system on a lattice, many 
details enter: lattice size and geometry, number of nearest 
neighbors (coordination number), interaction range, etc. To 
realize the square lattice LLV we adopt a typical Monte Carlo 
(MC) algorithm (details in [15,17]), namely (1) At every mi- 
croscopic step one lattice site is randomly chosen; (2) One 
of the nearest neighbors is also selected randomly; (3) If the 
original chosen site is Xi (X 2 ) and the selected neighbor is X 2 
(S) then the chosen site changes to X 2 (S) with probability k s 
(fci); If the original chosen site is S and the selected neighbor 
is Xi then the chosen site changes to Xi with probability k 2 ; 
Otherwise the system remains as it is; (4) Return to step 1. 

In the MC procedure the unit of time is chosen as 1/JV, 
where N is the total number of lattice sites (occupied and 
empty). For example, for square lattice, N = L 2 , where L is 
the linear size of the lattice. With this choice of micro-time, 
in one MC step all lattice sites are, on the average, scanned 
once. In Fig. lb typical behavior of the temporal evolution of 
the MC concentrations is shown. In particular, the concentra- 
tions of Xi is depicted on the full lattice of size L x L — 2 s x 2 s 
(solid line) and on a sublattice of size Ixi = 2 5 x 2 s (dotted 



line). Periodic boundary conditions are used in all simula- 
tions. It is clear that while on the sublattice the concen- 
trations show oscillatory behavior with added noise, on the 
entire lattice the oscillations shrink. Fig. 2a gives the typical 
evolution of a system starting from random initial conditions 
(Fig. 2a (t=0 MC)). As time increases the system develops 
local domains and each domain behaves as a local oscillator 
with specific characteristic frequency. Because the various 
domains have different phases, globally, no oscillation are ob- 
served, in contrast with the MF predictions [15]. Moreover, 
it has been shown [17] that the different species organize in 
local domains which present competing interactions and they 
have fractal boundaries. In this figure and hereafter the Xi 
particles are depicted in grey color, the X 2 in black and the 
empty sites in white. The initial condition was a homogeneous 
infinite lattice with equal concentrations of Xi , X 2 particles 
and empty sites S. Periodic boundary conditions are used. 
The fractal properties of the spatial structures can be used to 
measure the size of the local oscillators [17] and point out to 
a nonextensive formalism for the calculation of its entropy. 

To describe the temporal evolution of the entropy with re- 
spect to one of the species, e.g. Xi, we start from a given 
configuration, with specific initial conditions on lattice and let 
the system evolve according to the MC algorithm. The choice 
of the particular species does not play any role in the entropy 
calculations, since the model is cyclic and all the species are 
equivalent. As time increases the system passes through var- 
ious configurations which we record at regular temporal in- 
tervals. Let us call C(t) — {dj(t)},i — 1, ...L;j = 1, ...L the 
specific configuration of the lattice at time t, while Cij(t) de- 
notes the state of site at time t and 



Cij — 







if site is occupied by Xi 
if site is occupied by X 2 
if site is occupied by S. 



Within each configuration we introduce a set of M non- 
overlapping windows {Wi}, (i = 1, ...M), of size I x I which 
cover completely the lattice. The number of windows is 
M = n 2 = (L/l) 2 . Consequently C(t) = [jf =1 W>. Let us de- 
note with pi the probability that window i is occupied by par- 
ticles Xi. If ni(i, i) is the number of particles Xi inside win- 
dow i at time t and m (t) is the total number of Xi particles 
on the lattice, then pi (t) — ni(i, t) /m (t) . This probability set 
into Eq. (1) provides S q (t). For short times S q (t) scales as a 
nonlinear function of the time, while for large times depends 
non-linearly on the system size. This non-linear dependence 
on the system size is the basic indication of non-extensivity. 
The various values of q highlight characteristics on different 
length scales in the system. As an example, rare events are 
characterized by low values of pi. For q < 1, the term p? 
takes relatively large values and gives important contribution 
to the function S q . In contrast, if q > 1, then pi « pi and 
the contribution of rare events is negligible. 

It is well known that scaling behavior is proper to systems 
which present fractality. Especially in monofractals only one 
level of scaling is detected while in multifractal structures the 
different scales grow with different power laws. The S q en- 
tropy is then the appropriate measure of complexity because 
it addresses the complexity in different length scales by ap- 
propriate tuning of the q value. We study S q with different 
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initial conditions. Depending on the degree of organization 
of the initial state, the entropy may increase going to a more 
disordered state or decrease going to a more ordered state. 

Consider first the case of the "domain growth mode" , where 
initially the system contains particles X± , X 2 and empty sites 
S randomly distributed and with equal probability, as in Fig 
2a (t=0MC). The initial state is statistically uncorrelated. 
This is the state of maximum entropy. It is easy to calculate 
S q for this state as a function of the window size. Assume that 
the X\ particles are the information carrying sites, while the 
X2 and S are the medium. This assumption does not reduce 
the system complexity since all three species are equivalent 
for cyclic reactions. If the three species have the same con- 
centrations then the number of X\ particles will be on the 
average equal to xi(t = 0) = L 2 /3. If we divide the lattice in 
non-overlapping windows of size I x I then the average number 
of X\ particles in a box is I 2 /3 while the number of windows 
is n 2 = L 2 /l 2 . The initial values of p» are: 

» = I^ = (l) a CV0, (5) 

hence S q (0) = [{L/l) 21 - 1 ^^ — 1]/(1 — q). As time increases, the 
various species organize in domains with fractal boundaries 
and we call this mode "domain formation mode" . The entropy 
gradually decreases since the system evolves from a random 
to a more organized state: See Fig. 2a. In Fig. 2b the 
corresponding temporal evolution of S q is shown for different 
values of q. S q initially undergoes a few oscillations, while 
the system organizes in domains of homologous (identical) 
species, and then stabilizes into a lower entropy state. 

To calculate the entropy production rate at the nucleus 
growth mode, we start with a fully organized state consist- 
ing only of particles S and we include a nucleation droplet of 
infinitesimal radius r placed on the lattice. The droplet con- 
tains particles X\, X2 and S homogeneously and randomly 
distributed within the droplet area. As time increases the 
droplet grows forming spontaneously several rings of parti- 
cles X\, X2 and S sequentially. The widths of the rings 
shrink with the distance from the pure S region and when 
their width becomes zero the typical LLV fractal pattern ap- 
pears in the middle as can be seen in Fig. 3a. This type of 
spreading is called the "nucleus growth mode" because an ini- 
tially small droplet grows in size and finally covers the entire 
system. This 2-dimcnsiona (d — 2) growth leads to a reorga- 
nization of the species which at the beginning were randomly 
distributed within the infinitesimal droplet, while they even- 
tually present fractal patterns as in Fig. 2a. As time increases 
this typical pattern will cover the entire lattice and S q attain 
the values calculated in the previous case (see Fig. 2b). 

As seen in Fig. 3b, only the case q — 0.5 shows a linear 
increase with time during the entropy production duration; 
behavior is sub-linear for q > 0.5 and superlinear for q < 0.5. 
In Fig. 3c we see So. 5 for various system sizes. They all 
start linearly, coincide during the entropy production period, 
and saturate at different values, in accordance with Eq. (5). 
The q-value does not change with variations of the lattice size 
L, of the window size I, and of the initial concentrations of 
reactants within the original droplet. Note the difference in 
ordinate scale between Figs. 2b and 3b. If we inspect closer 
the L = 500 steady state of Fig. 3b, the entropy lines for 



t > 150 present fluctuations similar to the ones in Fig. 2b. 
An interesting data collapse is shown in Fig. 3d. In Fig. 4a,b 
the entropy of the LLV model, at the "stripe growth mode" 
is shown. This is a 1-dimensional growth (d = 1) mode. The 
initial state of the system consists of a stripe of randomly 
distributed X\, X2 and S embedded in a lattice containing 
only S particles otherwise (Fig. 4a). The entropy features at 
this mode are similar to the "nucleus growth mode" but now 
only q = produces the linear behavior, hence q depends on 
d. 

To explore the entropy production due to surface roughen- 
ing (or "roughening mode" ) , we investigate the case of an in- 
terface separating two stripes of identical particles. The setup 
of the system at the initial state is as follows: On a 500 x 500 
square lattice we create one stripe of size S x L — 50 x 500 
consisting only by X\ particles followed by a stripe of the 
same size but consisting only of X2 particles, while the rest of 
the lattice is covered by S, see Fig. 5a (t=0MC). In the Figs. 
5a, corresponding to times MC, 20 MC, 60 MC and 100 
MC respectively, the originally linear interfaces roughen and 
the stripes are deformed. The process is dynamical and all 
interfaces move to the right with the same average velocity. 
While the size of the stripes is on the average kept constant, 
fluctuations make them vary significantly. In fact, after suf- 
ficiently long times (depending on the width S and the size 
L of the stripes) the stripes will mix and the typical fractal 
patterns of Fig. 2 will reappear. The entropy increase due 
to the roughening of the interfaces is shown in Fig. 5b. Af- 
ter an initial increase due to roughening the entropy remains 
constant with statistical fluctuations around the steady state. 
This mode is clearly not adequate for extracting the physical 
value of q. 

In the current study the nonextensive entropic properties 
of the LLV model are examined. The special value of q which 
produces a linear increase of S q (t) depends on the dimension- 
ality d of the growth: q = 1 — 1/d (d — 1,2). These non-trivial 
values of q are in accordance with the appearance of fractal 
spatial structures observed in earlier studies for the same sys- 
tem; also, one expects q = 1 in the d — » 00 limit. Further 
studies (e.g., sensitivity to the initial conditions, multifractal 
function f(a), entropy relaxation, aging) are welcome. 
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FIG. 1. (a) The Mean Field Approximation concentrations 
xi(t) and x 2 (t) for (xi(0), x 2 (0), s(0)) = (0.5,0.4,0.1), and 
(ki,k2,k s ) = (0.9,0.3,0.1). (b) Monte-Carlo simulations; the 
solid (dotted) line corresponds to xi(t) over the full lattice 
(sub-lattice). 




FIG. 2. (a)Four different snapshots during the evolu- 
tion of LLV for random uniform initial conditions on a 
L — 500 square lattice at the "domain formation mode"; 
(fci,fc 2 ,fc s ) = (0.9.0.3,0.1). (h)S q (t) with I = 10. 



(t) 




FIG. 3. (a) Four different snapshots during the evolution of 
a system covered initially by S with one small mixed droplet. 
The system linear size is L — 500 while the droplet size is 
I = 8; (fci,fe,fc s ) = (1.0,1.0,1.0). (b)S g (t). (c)5 .b(*) for 
various lattice sizes, (d) Collapse of the Fig. 3c data. 
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FIG. 4. (a) Four snapshots of the LLV model at the "stripe 
growth" mode; (fci,fc 2 ,fc s ) = (1.0,1.0,1.0). (b)S (t) for 



500 x L lattices. 




FIG. 5. (a)Four different snapshots during the evolution 
of a L — 500 lattice containing initially two stripes of size 
50 x 500; (fci, k 2 , k a ) = (1.0, 1.0, 1.0). (b)S q (t). 
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